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Abstract 



We consider the evolution of perturbations to a flat FRW universe that 
arise from a "stiff source," such as a self-ordering cosmic field that forms 
in a global symmetry-breaking phase transition and evolves via the Kibble 
mechanism. Although the linear response of the normal matter to the source 
depends on the details of the source dynamics, we show that the higher- 
order non-linear perturbative equations reduce to a form identical to those of 
source- free Newtonian gravity in the small wavelength limit. Consequently, 
the resulting n-point correlation functions and their spectral counterparts 
will have a hierarchical contribution arising from this gravitational evolution 
(as in the source-free case) in addition to that possibly coming from non- 
Gaussian initial conditions. We apply this formalism to the 0{N) nonlinear 
sigma model at large and find that observable differences from the case 
of initially Gaussian perturbations and Newtonian gravity in the bispectrum 
and higher-order correlations are not expected on scales smaller than about 
100/i-^Mpc. 
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I. INTRODUCTION 



This paper attempts to combine several disparate streams of work in the study of cosmo- 
logical perturbations. Since the early work of Lifschitz [Q, Sachs and Wolfe and others, 
the evolution of large scale matter perturbations to a spatially homogeneous and isotropic 
expanding universe has been well-studied. More recently, this work has been refined to apply 
to a more complicated universe containing both radiation (with equation of state p/ p = 1/3) 
and non-relativistic matter {p = 0) constituents, using a gauge-invariant approach In 
addition, large A^-body simulations have been used to examine the fully nonlinear problem 
of gravitational evolution on scales much smaller than the Hubble radius, where the New- 
tonian limit of general relativity is sufficient . On the other hand, certain aspects of the 
full non-linear evolution equations do remain amenable to a more analytic approach. An 
area of work that has received considerable attention in the past decade is the examination 
of higher-order corrections to the linear Newtonian equations of motion for the matter in 
the universe within the framework of perturbation theory [[7|-[T0|]. These higher-order cor- 
rections are intimately related to the evolution of the multi-point correlations of the mass 
(and therefore, of the galaxies), such as the three-point function, ^3 cx (5(xi)5(x2)5(x3)), 
the skewness, or, in fourier-space, the bispectrum, B oc {6(ki)6(k2)S(ks)) ■ As increasingly 
large galaxy surveys become available, observational information on these moments of the 
galaxy distribution has been extended to scales sufficiently large that this perturbative ap- 
proach is expected to be reliable. In the case of an initially Gaussian density field for which 
all reduced higher-order correlations vanish, this "quasi-nonlinear" gravitational evolution 
leads to a scaling hierarchy of correlation functions: ^„ cx ^2"^; where is the volume- 
averaged n-point correlation function, and the constant of proportionality depends weakly 
on the power spectrum or two-point function of the perturbations. Even with an initially 
non-Gaussian distribution, it may be possible to distinguish the primordial component from 
that due to gravitational evolution on large scales. 

The standard cold dark matter (CDM) model for structure formation invokes a mecha- 
nism such as an early epoch of inflation to generate primordial adiabatic density fluc- 
tuation. Inflation results in a spatially flat universe overall; quantum fluctuation of the 
inflaton fleld lead to an initially scale-invariant spectrum of density perturbations (with 
details possibly depending upon the speciflc model of inflation) which evolve thereafter un- 
der the influence of gravity. In the simplest models of inflation, this initial distribution of 
perturbations is Gaussian, so the hierarchical results noted above obtain for the resulting 
correlation functions. The speciflc mechanism which drives inflation is inextricably linked to 
the fundamental particle physics model of the universe. Of course, in addition to providing 
a mechanism for the generation of perturbations, inflation also has the advantage of solving 
the horizon and flatness problems which are otherwise left unexplained. 

In another class of theories, topological defects (or some other classical fleld conflgu- 
ration) such as cosmic strings, domain walls or textures act as a continual source for the 
creation of density perturbations []TT|. These structures are the remnants of a symmetry- 



breaking phase transition of a cosmic fleld and are thus, like inflation, a cosmological relic 
of the underlying particle physics. In these scenarios, the fleld is initially laid down ran- 
domly before the phase transition, on scales larger than the Hubble radius at that time. 
Once the symmetry is broken, the fleld tries to align itslef in order to minimize its energy 
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density. However, this alignment can only occur coherently on scales where the field has 
come into causal contact with itself — within the Hubble volume. As the universe expands, 
the Hubble volume increases and the field orders itself on larger and larger scales — the Kib- 
ble mechanism |T^. If the symmetries of the field and its initially random configuration 
require it, topological defects may result. For example, a broken discrete symmetry like Zn 
produces domain walls, a broken 0(2) or f/(l) symmetry can result in (gauge or global) 
cosmic strings, 0(3) in monopoles, and global 0(4) in global texture. The breaking of a 
global O(A^) symmetry with > 4 does not lead to topological defects, but does result in 
spatial field gradients and consequently to perturbations to the energy density. Unlike the 
infiationary scenario, the Kibble mechanism quite generically produces density fiuctuations 
with an initially non-Gaussian distribution. (Of course, these theories do not provide a 
natural solution to the horizon and fiatness problems and in this case it is usually assumed 
that the universe begins in a perfectly homogeneous and isotropic state with f2 = 1 before 
the symmetry-breaking phase transition. Thus, in these models one relegates the solution 
of these puzzles to initial conditions, or to an earlier infiationary epoch driven by a field too 
weakly-coupled to produce the density perturbations responsible for large-scale structure.) 

Within the class of scenarios that create structure via the Kibble mechanism, there is 
a further conceptual split between global and local (or gauge) symmetries. In a local the- 
ory, the large-scale gradient energy of the field is compensated by the gauge field so all the 
field energy is concentrated in localized defects. Thus, only true topologically stabilized 
configurations can result. In global theories, on the other hand, non-topological configura- 
tions (textures) can result that are nonetheless long-lived because they require energy to 
"unwind" a configuration by forcing it off its vacuum manifold. In large-A^ models, these 
field configurations persist simply because causality constrains them to align only on scales 
smaller than the Hubble radius, so field gradients persist for approximately a Hubble time. 
(Also, a new class known as "semi-local" defects has been studied, in which a gauge theory 
admits defects which are stabilized by the dynamics of theory, not the underlying topology 
of the symmetry groups ||13| .) 

In this paper, we shall modify and extend aspects of a formalism that has been developed 
by Veeraraghavan and Stebbins to study the perturbations due to a "stiff source" such as 
these cosmic field configurations formed by the Kibble mechanism. A "stiff source" evolves 
in the homogeneous and isotropic background metric of the universe; the back-reaction of 
the metric perturbations onto the source is considered to be negligible. We shall, however, 
explicitly account for compensation: the initial response of the matter fields to the stress- 
energy of the stiff source. We shall extend previous work to allow the perturbations of the 
matter and radiation fluids in the universe to enter the quasi-nonlinear regime and examine 
the modifications to the perturbative solution to the equations of motion that result. This, 
in turn, allows us to study the resulting correlation functions and may modify the scaling 
hierarchy in such scenarios. 

Using this formalism, we shall concentrate on a specific class of models in which the 
Kibble mechanism is responsible for the initial generation of density perturbations, the 
nonlinear 0{N) sigma model ||T5|-[T7|. These models arise as the low-energy limit of the 
breaking of a global 0{N) symmetry to 0(A^ ~ !)• Although the analytic calculations we 
shall perform rely on the large- A^ limit, in which there are only spatial gradients, but no 
topological defects, we can also extract some information about the behavior of this theory 
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for small where defects play a role. 

The plan of this paper is as follows. In Section |T|, we develop the equations of motion for 
the stiff source in the background metric and of the matter fields in the perturbed universe, 
using the longitudinal gauge, and develop a perturbative expansion about homogeneity and 
isotropy. In Section |T|, we apply this formalism to the nonlinear sigma model, and examine 
the higher-order correlations of the matter. Finally, we present our conclusions. In an 
appendix, we show some useful results for the distribution of the fields in an 0{N) model. 



II. COMPENSATED PERTURBATIONS 

A. Perturbations in the Longitudinal Gauge 

We will consider linear metric fluctuations about a homogeneous FRW background. We 
work in comoving conformal coordinates with metric signature (+,—,—,—) and assume a 
spatially flat background metric of 

ds"^ = a^rjuiydx^dx'^ = {^drf — Sijdx^dx^^ , (1) 

where 77^,^ is the Minkowski metric. Here, t] is conformal time, related to proper time 
by dt = adr]. Throughout, we shall use a prime to denote the derivative with respect 
to conformal time, and define a conformal expansion rate 7i = a' /a. The unperturbed 
Einstein equations in the fiat (f2 = 1) Freedman- Robert son- Walker (FRW) universe with 
mean background density p and pressure p and vanishing cosmological constant (A = 0) are 

H = -—a p; — = -—a [p - 3p). (2) 
6 a 6 

In a matter-dominated universe {p = 0), the scale factor a ccrf and in a radiation-dominated 
universe {p = p/3), a cc rj. The Hubble constant is Hq = lOO/ikm sec~^ Mpc~^ and the 
present Hubble radius is H^^ = 3000 h~^Mpc. We normalize the conformal time by setting 
?7o = 2-ff(7^ = 6000 h~^Mpc today. We will write the perturbed metric as 

Qfj^v = a^iVi^v + ^^,1,). (3) 

In the longitudinal gauge (/loi = 0) to first order in /i, g'^^ = a^'^{rj^y—hfj,y) or g^^ = a^(l — /iqo) 
and g''^ = —a~'^{6ij + hij). 

In the usual longitudinal (often called conformal-Newtonian) gauge analysis, only scalar 
perturbations are considered, in which case the metric perturbations are determined by two 
scalar variables ("potentials"), /iqo = 20 and h = 6ijhij = Gip- In this paper, we will in addi- 
tion allow nonzero vector and tensor perturbations, for which h^j = h^j — 6ijh/3 is nonzero. 
Most analyses of density perturbations in an expanding universe have used the synchronous 
gauge (/loo = = 0) [0 (although a gauge-invariant approach has recently become popular 
0-^). We have chosen the longitudinal gauge in order to more easily compare our results 
with the perturbative Newtonian equations of motion (and thereby gain insight into higher 
order correlations which have previously been studied only in the Newtonian limit); the 
price of this is retaining a nonzero component /ioo- More importantly, in the longitudinal 
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gauge we can consider situations in which the density perturbation amphtude S = Sp/p is 
large, while the metric perturbations are still small. 

A few more words about our approximation scheme are in order. Because we shall only 
compute the Einstein tensor to first order in the metric perturbation, the perturbed Einstein 
equations 5(7^,^ = 8nG6Tfj_i^ (where here 6T^u represents the perturbation to all the stress- 
energy in the universe, including matter, radiation and any sources such as those discussed 
below) can only be considered first order — in the metric perturbation. Thus, we must 
ignore terms like h^^, ■ dT^p for consistency. However, this does not determine the form 
that the stress-energy must take: it may contain v"^ or v5 terms. (Below, we shall calculate 
the four- velocity to order f^.) Thus, these equations can be valid for "large" values of the 
density perturbation amplitude 5 ~ 1 as long as we still have f ^, h^^ -C 1. Note that 
this is a gauge- dependent statement. In fact, performing the same exercise in a comoving 
synchronous gauge, we find that h' = 26' to all orders in the matter perturbation variables 
for pressureless, nonrelativistic matter, so this scheme would not be successful — we can only 
use synchronous coordinates until streamlines of the matter flow intersect and caustics form. 
The advantage of the synchronous formalism is that one can define a comoving gauge where 
the nonrelativistic matter component has zero velocity for all time {i.e., as long as the 
coordinate construction is consistent). In the longitudinal gauge we are considering here we 
still have to calculate velocities. However, the gravitational potential (metric perturbation) 
in this gauge is in general suppressed compared to the density perturbation by the square of 
the size of the perturbation relative to the Hubble radius, as in the usual Newtonian Poisson 
equation: ~ ~ [X/H~^)'^S for perturbations of scale A. (Of course, this equation is 
itself only applicable for scales smaller than the Hubble radius, beyond which relativistic 
corrections become important as in e.g., Eq. ([T8| ) below.) 

On scales smaller than the Hubble radius, the metric perturbation remains small even 
when the matter perturbation becomes nonlinear. In particular, we expect this approxima- 
tion to be valid as long as h"^ <^ 5; with first order quantities and the Poisson equation to 
relate the potential to the density perturbation, this is equivalent to {X/H~^)^ -C 1/5. As 
expected, our approximation will continue to hold on scales much smaller than the horizon 
as long as 5 ^ 1. Conversely, for sufficiently small scales the linear metric perturbation 
approximation holds for even larger values of the density perturbation. As long as the av- 
erage source stress-energy is negligible compared to the fluid components and as long as we 
continue to ignore backreaction, this approximation holds in the same regimes as the normal 
Newtonian limit. 

In the longitudinal gauge, hij does not contain a scalar part as it does in the synchronous 
gauge, so in addition to being manifestly traceless, it satisfles didjhij = 0. Although the 
matter variables become more complicated when we allow higher-order terms to enter the 
equations, the metric perturbation is still usefully decomposed into geometrically distinct 
parts {i.e., scalar, vector and tensor) Q]. Normally, only the scalar part of the peculiar 
velocity, which can be expressed as the gradient of some velocity potential, appears in the 
equations of motion for other scalar quantities like the density perturbation S. With higher 
order terms, however, products like f *5 or v^v^ decompose as a vector and symmetric tensor, 
respectively, with scalar, vector, and (for v'^v^) tensor parts that are not simply related to 
the scalar and vector parts of the original f*. Thus, for example, 6, although a scalar, can 
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appear in higher-order equations containing the vector perturbations to the metric, and the 
vector part of can occur in the scalar equations of motion {e.g., Eq. (|9al) which contains 



As in (who instead calculated in synchronous gauge), we assume a universe filled 
with a multi-component perfect-fluid stress-energy T^^, = Y,nT(l^y, where n labels the fluid 
component, as well as a "stiff source" with stress-energy Q^p. The source makes a negligible 
contribution to the spatially-averaged stress-energy and is covariantly conserved with respect 
to the background metric a^rj^^: 

e'oo + n (Goo + 6) = 5^60.; (a^Bo,)' = a^d^Q,,. (4) 

Here 6 = Qu is the spatial trace of the source stress energy. 

Again, we do not consider the back-reaction on the stiff source of the metric perturba- 
tions. In principle, we could calculate it and its effects in an iterative, perturbative fashion, 
writing the full source stress-energy as a sum of a part propagating in the FRW background 
as in Eq. and another part propagating in the full perturbed metric. This perturbation 
would, in turn, be fed back into the field equations for the metric perturbation. This would 
result in new terms of order {d^hap)Q-ys, still at linear order in the metric. Unless either 
the source stresses or the metric perturbation are large, these terms should be negligible in 
comparison to those involving the linear fluid perturbations. 

The bulk of the universe is filled with several perfect fluids with individual stress-energies 

T(ny = (Pn + ^Pn + Pn + Spn)KUnu - {Pn + ^nWu, (5) 

where n labels the fluid component, pn and pn refer to the background density and pressure 
of each fluid, and 5p„, 5pn to their respective perturbations. The four- velocity of each fluid is 
= vy + , normalized to u^u^ = +1 with this signature. In the expanding background, 
the fluid components are at rest, with four- velocity (suppressing the fluid label n) 

= {a~\ 0, 0, 0), Up = (a, 0, 0, 0) (6) 

and the velocity perturbations are 

Su^ = ^«~^(^^ - ^oo), Suq = ^a{v'^ + hoo), 



6u' = -a-^6ui = a~\' (^1 + ^v^^ (7) 

This defines the peculiar velocity 3- vector, f* = u'' /vP. We have kept terms of order f^, and 
ignored terms of order v ■ h, for the consistency of our approximation scheme. (Note that 
the form of the f ^ terms is exactly as one would expect from the usual special-relativistic 
four velocity expanded to this order.) 
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B. Equations of Motion for the Matter and Metric Fields 



We will assume that each of the perfect fluid components is separately conserved. Except 
when there is significant energy transfer between the components {e.g., at the epoch of 
recombination) this is an excellent approximation. To first order in the metric perturbation, 
the equations of motion for a single fluid component, V^jT^^y ~ 0; become 



5' + v^do [5(1 + C)] + (1 + u» + 5(1 + 0) + V • [(1 + w + 5(1 + C)) v] 



+ m 



v'^{l + w + 5(1 + C)) {w + -) + 5(C -w)+ v'^{l + w){w- c^) 



(8a) 



and 



Hv' [{1 + W + 5(1 + 0) (1 - M + 3(1 + w){w- c')] + [5(1 + C)]' v' + 

(1 + w + 5(1 + 0) v" + dj [{1 + W + 5(1 + 0) v'v^] + 9,(C5) = -(1 + w)d4 (8b) 

where 5 = 5p/p, C = ^p/^P, the fluid sound speed = dp/ dp, w = p/p, and we have 
suppressed the subscript n on all fluid quantities. To linear order, (n = and we shall often 
take Wn = Cn = as well. In the background, the fluid density evolves according to the 
zeroth-order equations of motion, p'^ = —3(1 + Wn)pn- Notice that these equations do not 
depend exphcitly upon the stiff source. 

The perturbed Einstein equations 5G^y = SnGiJ^in) ST(^n)fiu + ©m^^) become 

= m (^' + n<P) + ATiOa^ Pn K + vl{l + Wn + 6n{l + Cn))] + 47rGeoo (9a) 



-nW(t) - AT^Ga^ Pn^n {l+Wn + 5„(1 + C„)) + 47rGe 



Oi 



(9b) 



A-KGTi. 



^" + -V2(0-^) 



+ ni2f + (f)') 



+ AnGa^ ^ p„ ISnCJij + (1 + w„ + 5„(1 + Cn)) vlvi] + ATrGOi 



(9c) 



where 



(10) 



is the trace-free part of the spatial Einstein tensor in this gauge. We have written the 
equations such that the form of the stress-energy pseudotensor Tap is manifest. By Eqs. (H), 
the pseudotensor manifestly obeys 
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d^T^" = 0; 



(11) 



In effect, the pseudotensor can be written as 



grav 



(12) 



where T^^ and Q^iy are the ffuid and source stress-energies, respectively, and the remaining 
term, t^^^ represents the "stress-energy of the gravitational field." As we see in Eqs. (H), 
this term contains those parts of the Einstein tensor that are either nonlinear in the metric 
perturbation (although these terms do not appear in our approximation) or those that would 
vanish for a constant scale factor a' = 0. 

The pseudotensor is conserved with respect to a Minkowski background; it thereby em- 
bodies the concept of a conserved stress-energy for the entire system including the "gravita- 
tional energy" as well as the matter and source stress-energy |l^Jl^JT^ , the combined system 
propagating on a flat, Minkowski background. However, no individual fluid or gravitational 
components of the pseudotensor can be singled out as being separately conserved in this 
Minkowski space; moreover, its definition is gauge- dependent. It is useful for setting up 
the conditions of compensation on scales larger than the Hubble radius: it can describe the 
flow of energy between matter and "gravity" which only occurs causally, on scales within 
the Hubble volume. Note also that to this order, the metric perturbation only appears in 
terms containing factors of Ti. — in a true Minkowski background [i.e., not expanding), the 
pseudotensor only differs from the "real" stress-energy tensor at second and higher orders 



in the metric perturbation ||T8|,|T9|. 

The space-space component of the Einstein equations may be rewritten as a trace and 
traceless part: 



2^- 

a 



didj 



+ n{2ij' + 0') + ij" + -V2(0 - ij) 



AtiGo? ^ p„ 



SnCn + -{l+Wn + 5„(1 + Cn)) 



+ 



4nG 



- V^) + ^ [v^/iij - djdkhik - didkhjk - h'-j - 2nh[^ 

= StcGo^ ^ P„ (1 + M;„ + 6n{l + Cn)) - \SijV, 



e 



(13a) 



STiGBij. (13b) 



Here, Eq. (|13a|) is the spatial trace of Eq. (|9^) (AirGrkk) and Eq. (|13b| ) comes from subtracting 



off that trace from Eq. (p^ ) (rij — SijTkk/'i)- Also, Qij = Qij —6ijQ/3 is the spatial trace-free 
part of the source stress energy. 

(Of course, the conservation Eqs. (||) are not independent of the Einstein equations, 
which imply that the sum of the individual fluid stress energies is covariantly conserved 
via the Bianchi identity. In an approximation that we will consider several times later, 
a universe with only one fluid in addition to the stiff source, the conservation equations 
are redundant and can be derived from the Einstein equations by using the conservation 
formulas for the pseudostress-energy (Eq. (pT])) and eliminating the source terms Q^u by 
using their background equations of motion, Eqs. (^.) 
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Obviously, these equations (P, H, 0) cannot be solved analytically for a general equation 
of state, with general initial data and an arbitrary stiff source. However, we can explore 
various limits and approximations that are useful in different theories of structure formation. 

C. Correspondence with Newtonian Equations of Motion 

If we consider a universe filled only with pressureless matter {w = = ( = 0) and no stiff 
source {Q^u = 0) we have the conditions of the standard matter-dominated universe. If we 
further assume the limit of scales small compared to the Hubble radius and small velocities, 
f <C c, then we have the situation usually described by the Newtonian limit. In the equations 
of motion for the matter (Eqs. (H)) and the Einstein Eqs. (^|13]), we drop terms like d^v"^ in 
this limit but we keep terms like div"^ because we assume that the matter component may 
vary over small scales. Note that we do not assume that the density perturbation 6 is small, 
so we do not ignore terms that contain 6 in spatial derivatives, even when multiplied by 
small quantities like Ti.^ or v. In this case, the spatial Einstein Eqs. (|13D only describe the 
evolution of the scalar components of the metric, and ip; the tensor-mode equations drop 
out as expected as these components evolve completely independently. The i ^ j Einstein 
Eq. ( |13b| ) then requires cj) = ip. The trace of the i — j Eq. ( |13aD gives 

a" \ 

2 n^\(p + met)' + 0" = (14) 

a J 

with solution (j) = const (plus a decaying mode proprotional to r]^^). With these results, 
the — Einstein Eq. (|9a|) and the stress-energy conservation Eqs. (|^) (which are not 
independent from the Einstein equations) become 

^n^S = V20, (15a) 

5' + V ■ [(1 + 5)v] = 0, (15b) 
T^v + v' + (v ■ V) V = -V0, (15c) 

where we have used the second equation to eliminate 6' in the final equation. 

These equations are just those that come from a purely Newtonian analysis of density 
perturbations in the fluid limit: the Poisson equation, the continuity equation, and the 
Euler equation. It should not be surprising that this approximation scheme has yielded 
these results: in Newtonian theory, the potential ~ f ^, so it is suppressed to second order 
compared to the matter variables when the velocities are small, and ignoring terms like 
(f)6 and (pv should be sufficient to this level of approximation. This approximation can be 
made more precise by expanding separately in two small parameters: the size of the metric 
perturbation h and the ratio of the length scale of the perturbation in question to the Hubble 
length. This expansion generates the Newtonian and post-Newtonian approximations in a 
cosmological setting |[18| , |20[| . 

If we consider a simple perturbation expansion in both the matter velocity v and density 
perturbations 6 (with both quantities the same order), and an initially Gaussian density field 



9 



{e.g., from inflation), these equations have been shown to produce a scahng hierarchy of cor- 
relation functions in the quasi- nonhnear regime [0,||,|lO|: oc ^2~^5 where ^„ is the volume- 
averaged n-point correlation function or moment of the distribution {e.g., ^3 is the skewness). 
For the unaveraged multi-point correlation functions ^„(xi, . . . , x„) or their fourier-space 
spectral counterparts Pn(ki, • • • , kn), the n-point functions are proportional to sums of ap- 
propriate symmetric products of {n — 1) two-point functions. For the simple case of the bis- 
pectrum, S123 = Q(-Pi-P2 + -Pi-P3 + -f2-P3) as in Eq. ( pD| ) below; the trispectrum is a product of 
three power spectra, T1234 = Ra[P{ki)P{\ki+k2\)P{h)+sjm.]+Rh[P{h)P{k2)P{h)+sjm.], 
where "sym." represents the further terms of the same form necessary to keep T symmetric 
in its arguments. Here, Q, Ra, Rb and their generalizations to higher n-point functions 
are constants that depend on the configuration of the kj and possibly on the initial power 
spectrum. 

For this hierarchical result to apply, two additional assumptions are necessary: the de- 
caying mode component of the solution for the density perturbation must become negligible, 
and the velocity-field must be curl-free, so it can be expressed as the gradient of some po- 
tential; any leftover vortical component of the velocity decays away proportional to a~^, so 
this assumption is justified at late times. 

It is interesting to examine how these equations are modified even in a universe with only 
a matter component when relativistic effects are included. That is, we no longer make the 
Newtonian approximation of small velocities and length scales that was used to derive Eqs. 



(p!5|). Even to linear order, these equations have corrections for perturbation wavelengths 
larger than the Hubble radius (A ^ H~^, which we shall hereafter refer to as superhorizon 
scales). The Poisson equation simply becomes the non- relativistic matter {w = 0) version 
of Eq. ( pal) with (p = ip, and the continuity equation gains a potential term (relaxing the 
assumption (p = const): 

^n^s = v^(/) - 3n(t>' - m^cp (i6a) 

30' = 5' + V ■ [(1 + (5)v] (16b) 
while the Euler equation is unchanged. The linear solution to these equations is well known. 



and is presented in the following section, Eqs. (|T8]). The full nonlinear equations also have 
corrections of order f^. However, we expect that when these terms are significant, our 
approximation that the metric perturbations are still of linear order will break down. 



D. Full Equations in the Linear Regime 

One useful application of Eqs. (|8|-|l3|) is to a universe consisting only of pressureless mat- 
ter, radiation and, possibly, a stiff source, in the linear regime (where now we are considering 
perturbations linear in everything: 5, and h^y). Because these equations are linear, the 
full power of previously developed techniques can be applied. In particular, we can use the 
geometrical decomposition of these equations to study the scalar, vector, and tensor modes 
separately. Density fluctuations in particular only couple to the scalar mode of the metric 
perturbation to this order, (p and ip, as well as the scalar parts of the peculiar velocity and, 
if present, the stiff source. In a universe with non-relativistic matter {w = = C, = and 



10 



radiation 
become 



w 



( = 1/3) fluids, the linear Einstein Eqs. ([9a|j9B| , [l3D for tlie scalar modes 



00 



i^fj' + H(j)) - ATtGdiQoi 
1. 



-^) 



3 ^ 



-SttG 



3 ^ 



(17a) 
(17b) 
(17c) 
(17d) 



where f2m and fi^, respectively, represent the ratio of the matter and radiation densities to 
the (time-dependent) critical density, with fi^ + fir = 1 by assumption. Q^j^ is the scalar, 
trace- free part of Qij, and can be written as above in terms of a second-order linear operator 
acting on a unique scalar function s{7],x.), determined, for example, by expanding Qij in 
terms of appropriate harmonic functions. (In particular, Qij ^ SijQ/3 on superhorizon 
scales, where the differential operator that projects out the s component approximately 
vanishes, as do all spatial derivatives.) Thus, the final equation tells us that (f) — il) = —SnGs 
up to a spatial constant (and in a universe with no stiff source, = as usual). Even without 
a stiff source, these equations cannot be solved exactly in the presence of both matter and 
radiation fluids except in various limits. In a matter-dominated universe {Qm = 1), the final 
two equations give (p and ip in terms of the source, and the first two in turn give the density 
and velocity. 

Because we are again considering linear perturbation theory, the complete solution in the 
matter-dominated era is just given by the sum of the usual pure matter solution (vanishing 
Qr and Qfiu) and a "particular" solution generated by the Qf^^ terms. The source-free matter 
solution is: 



1 
6 



= Ci(x) + C2(x)r7 ^ ^ const 
r^V^Ci - 12Ci + UV^C2 + I8C2) V 



-5 



6 



120. 



;i8a) 



;i8b) 



Superhorizon scales can be defined by the condition that r^V 0; in fourier space where 
V — >■ ~ 1/A, this becomes the appropriate ki] 0, since we also have t] ~ H~^. On 
these scales, 6 ~ —20 const, so perturbations do not grow, while on small scales the 
usual Poisson equation applies. In this case, the equations can be rewritten as the linearized 
version of the Newtonian case considered above in Section |II C| . 

In the presence of a source term, we must add the appropriate particular solution to this 
homogeneous solution. After setting — -0 = —StcGs as discussed above, the solution for 
the potential is 



+ 1 [ XVdv- [ rj^X dm xiv, x) = SvrG (lo - 2Hs' - s" + ^ V^s 



(19) 



5 7 ' 5 ' 7 ' " / 3 

Given these solutions for the scalar potentials, the density fluctuation can be computed 
immediately. In the presence of a source term, the question of initial conditions for these 
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fluctuations is crucial. We wish to express the fact that, at the time of "source creation" 
(e.g., a phase transition) the universe is initially homogeneous, and the matter variables can 
only respond causally to the source stress-energy — that is, the initial anisotropics of the stiff 
source stresses must be compensated by the matter fields. Within the Hubble volume, this 
produces actual perturbations in the fluid component; over superhorizon scales, the matter 
fields cannot vary coherently due to causality, so the universe must have zero curvature 
outside the Hubble volume. This isocurvature nature of these perturbations is embodied in 
the Minkowski-space conservation of the stress-energy pseudotensor t^'^ whose components 
were given above. In the initially homogeneous universe, the pseudoenergy Too = = as 
well as = 0. As is manifest from the Einstein Eqs. (j^), the evolution equations for the 
pseudotensor are given by 

^00 = diToi] = djT.j. (20) 

When the phase transition, a Poisson random process, occurs, we expect both of these 
quantities 0^^/ and r^j^ to be initially uncorrelated with themselves and thus gain white noise 
power spectra on superhorizon scales. That is, we Fourier decompose all quantities as 

/(k)= |rf3x/(x)e^'^-^ (21) 

and we shall write / ~ A;"; this means that the power spectrum l/^p oc A;^". Thus, we expect 
~ and r^^, ~ k^. The i — j components only occur in spatial derivatives, so these 
white noise spectra should obtain for these components for all time as long as fcr] ^ 1 {i.e., 
while the modes are outside the Hubble radius and out of causal contact with themselves), 
where we are considering the contribution from logarithmic intervals about a wavenumber 
k. From the evolution equations, we then see that r^i ~ k and Too ~ fc^, whereas 9oi ~ k 
and 6oo ~ k^ (since the Gqo equation is dominated by the H term). Thus, the quantities 
Too, TQi, and Goi should remain negligible for ki] <^1. A similar argument using the Einstein 
equations reveals that the density and metric fluctuations should initially have white noise 
spectra on superhorizon scales, but that velocities should fall as v ~ fc. 

Therefore, the superhorizon scale density perturbation is simply given from the — 
Einstein Eq. (^) by the relation Too = 0: 

^ nj^ = -\n-^eoo - 2H- V - 20 (22) 

n "J 

Here, this equation holds for a mult i- component universe; below, we shall specialize to a 
universe dominated by a single component with Pnl Pn = and n„ ^ 1. On large scales 



ip' ~ (5^/3(1 + Wn) (from Eq. (|8a|)) and ~ const. For the primary fluid component, this 



gives the superhorizon evolution equation for 5, 



^ + o^nl = -|^~'0oo - 20 (23) 

3n(l + Wn) 3 



For a universe with scale factor a oc rj' 

1 + w 



S{r], x) = ±^r7-3"(i+-)/2 / dr] eoo(r/, ^)r]^»(^+-^y^+^ - 20(r/, x) (24) 

a J 
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In a radiation-dominated universe (w = 1/3) this gives 5r- (In order to calculate the super- 
horizon matter perturbation in this case, we define the entropy perturbation a = 3Sr/4: — Sm- 
If as causality requires there is no initial entropy perturbation on superhorizon scales, 
6m = 3(5r/4 and it is unnecessary to explicitly worry about the details of the several compo- 
nents.) 

In particular, if the stiff source obeys a scaling relation 0oo oc r/"^ (see Section PD , then 



5 = -7T^r/2eoo-20. (25) 
Both of these terms are time-independent, so the perturbations do not grow with time outside 



the Hubble radius. The same result was derived by Davis et al. in the synchronous gauge [jT6 
This leads to the usual scale-invariant Harrison-Zel'dovich power spectrum: 5p/ p ~ const at 
horizon-crossing. In the nonlinear sigma model to be discussed later, the matter-dominated 
era evolution does obey 6oo oc r^"^, but in the radiation era Boo oc ln(?7/r/i), where 
Tji is the conformal time of the symmetry-breaking phase transition, so the evolution of 
the density perturbation is modified by a (divergent) logarithmic term. Thus, the initially 
white-noise spectrum (constant amplitude on all scales at one time) has been transformed 
into a scale-invariant spectrum (constant amplitude at horizon-crossing) by the gravitational 
action of the stiff source. 

Compensation thus insures that the pseudoenergy £ vanishes on superhorizon scales for 
all time. This fact in turn gives the initial condition for the perturbation on a given scale 
at Horizon crossing. (Note that this is not the case for the usual primordial adiabatic per- 
turbations discussed in CDM models, where the initial perturbation spectrum is considered 
as a given before the equations are integrated.) 

In order to calculate the evolution of density fluctuations well inside the horizon, it is 
easiest to use the first of Eqs. (0), Tqq = diTQi and the Einstein equations which define the 
components of r^,y, Eqs. (^. Here, it is assumed that the universe contains both radiation 
and matter fluids. After some algebra to eliminate the matter velocity and the potential ip, 
this component of Eq. (pO]) reduces to 



ml + 7i^6l - '^n'nmSm - m\i - fi^)5, + m\am - 2)</) + 37^V 

= A'kG (S.Goi - e^jo) (26) 

With the initial condition £^ = due to compensation, in a matter-dominated universe, this 
becomes simply 

^H^(5 + 2(1)) + m' = AnG p dr] diOoi (27) 

where we have ignored the contribution of Oqo to the total energy density inside the horizon. 
The solution is 

5= ^//-'y rfr/ir/^y drj2dieoi{v2)-Qv''j dr]r]''(l) + Kr]-^ (28) 

where we have explicitly included a decaying mode oc r]^^. A similar equation was derived in 
the synchronous gauge by Davis et al. |T^ . In the synchronous gauge, the term involving is 
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not present and the equation can be solved directly for the density perturbation. This term, 
however, is exactly as one would expect in performing the gauge transformation from syn- 
chronous gauge to longitudinal gauge for a scalar quantity 5p\ 6i = 6s—3H{l + w)a~^ J drj acf) 
1^, where / and s refer to quantities in the longitudinal and synchronous gauge, respectively. 
In the longitudinal gauge, unfortunately, things are more complicated and we must use our 
solution for the potential from above subject to the initial condition of S{kr] ^ 0) = 
for a complete solution. However, on extreme subhorizon scales, we can in general neglect 
terms like 7i^0 in the above equations — the longitudinal and synchronous coordinates nearly 
coincide. For example, in taking the Newtonian limit of the ^Goo Einstein equation in a 
pure CDM universe, we drop such a term in order to derive its nonrelativistic limit, the 
Poisson equation. Thus, the growing mode solution reduces to solely the first term above. 
In this case, and also assuming that the source has the power law form / dt] diOoi oc rj'^ (as 



in Section |T| below, the growing mode solution is 

47rG 



2(5 + 7) 



dr]diQoi- (29) 



Davis et al. derived this result for the specific case of 7 = in the synchronous gauge []16 
in which case the solution is of the same form as the normal pure-CDM growing mode 
6 oc a (X 7]"^ (which would be the solution here in the case of a nonzero initial value of the 
pseudoenergy: 6 oc S^rf.) 

In a radiation-dominated universe, the situation is more complicated. Note that we are 
interested in the evolution of matter perturbations in this situation. Thus, we take VL^ ~ 
in Eq. (pO) above, giving 



C + m'^ ^ A7iGn-\d,e^, - e[,o) (so) 

where the terms involving (f) and 5r have been ignored well inside the horizon. The solution 
is 

d7]2V2'j d7]iv!idiQoi-<doo)m^ (31) 

which include the usual constant and logarithmic terms, in addition to one due to the source 
evolution. Again considering the case of / drjdiQoi oc rf , the density perturbation due to 
the stiff source becomes 

8 = 4vrG'^^-^r/2 j dr, dSoi- (32) 

In the particular case that / drjdiQoi ~ const, (7 ~ 0) the response to the stiff source 
exactly mimics the usual isocurvature perturbation scenario: perturbations only grow inside 
the horizon and during matter-domination, 5 oc a (or logarithmically during radiation- 
domination) . If the source varies more rapidly, then corrections to this behavior may become 
important. 
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E. Higher Order Equations With a Stiff Source 



On superhorizon scales, the linear equations considered in the last section cannot be 
extended to higher orders in the matter variables within the approximation of linear metric 
fluctuations. Outside the horizon, large density fluctuation create large fluctuations in the 
metric, and we would at least need to treat higher orders in h^^i, as well as the matter 
variables (if not solve the full equations numerically). Well inside the horizon, however, 
the situation is different and the density fluctuation amplitude may be large for a small 
metric perturbation in the longitudinal gauge. Thus, we may still approximate the initial 
conditions as a vanishing pseudoenergy, but retain the nonlinear terms on small scales. 

Further, on small scales, we may again assume that the usual conditions of the Newtonian 
limit apply: potentials and velocities are small and slowly- varying. Of course, because of 
the presence of the source terms, it is important to at least check that these conditions 
still hold. From the linear solution of Eq. (P7D, we expect the potential to vary rapidly 
in regions of spacetime where the source (speciflcally, its scalar parts G and s) is rapidly- 
varying itself. For the sigma-model discussed below, the time scale of variation is always 
the Hubble time, so this assumption is justifled in that case. For models with topological 
defects, the regions around the defects themselves would generally have large stresses, and 
we expect this analysis to fail there. 

In order to compare with the usual Newtonian equations, the relevant equations are again 
the covariant conservation equations and the — Einstein equation, now supplemented by 
the traceless i — j equation to take into account the anisotropic stresses of the stiff stress- 
energy. On subhorizon scales, these give 

^n^S = VV, (33a) 

5' + V • [(1 + 5)v] = 0, (33b) 
+ v' + ( V - V) V = -V0, (33c) 
- ^ = -SnGs. (33d) 

(A rapidly- varying potential would add a term 3ip' to the right hand side of Eq. ( p3b|) .) The 
linear solution to these equations for stiff sources originating in a smooth universe is just that 
derived above. Although these equations do not explicitly contain the term Goi that appears 
in the above solution (Eq. (pQ])), recall that the source necessarily obeys conservation with 
respect to the background metric, so the two forms cannot be independent; they are related 
by the equations of motion for the source, Eqs. (^). 

Except for the presence of the s term, these are exactly the same equations as in the 
Newtonian matter-dominated case without sources, discussed above, Eqs. (ITH). Note also 



that s occurs in a "pure source" term — it is a function that is supphed from outside of this 
set of equations and the set of equations are linear in that function. 

This permits a particularly simple derivation of the higher-order equations of motion — 
i.e., the deviations from the linear solution above. We can write the various quantities 
as 

S = ^'^S + ^'^S + --- (34) 
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where the superscript represents an n-th order quantity {i.e., 0[{^^^S)"']) . Once the hnear 
equations have been subtracted off, the resulting equations give differential equations for the 
higher-order terms in terms of lower-order ones. These equations will have exactly the same 
structure as in the pure Newtonian case (Section |II Q above); the only dependence upon 
the stiff source will be implicitly through the ^^^6 solution. That is, we can write the linear 
solution as a functional of the stiff source, '^'5 = As usual in perturbation theory, the 

higher-order quantities only explicitly depend on the lower order '"'5: 

= [WS, . . . , ("-^'5] , (35) 

where the dependence on the source is implicit in '^^5, and thus in all higher-order terms as 
well. 

In the Newtonian case, it is usually assumed that the initial (linear) fluctuations have 
a Gaussian distribution; in the presence of a stiff source the distribution of the density 
fluctuations depend crucially upon the distribution of the source through the linear solution. 
However, the distribution of a stiff source is not a general property, but depends upon the 
specific model. As noted above, if / drjdiQoi ~ const, the matter perturbation behaves just 
as it does under the usual Newtonian evolution. Then, the problem is equivalent to the 
Newtonian one with the additional possibility of non-Gaussian initial conditions [BlI . 



III. THE NONLINEAR SIGMA MODEL 
A. Background Evolution 

Although the derivation thus far has been completely general and valid for any stiff source 
that obeys conservation with respect to the background (except as noted), most sources do 
not have a simple analytic form that can be "plugged in" to these equations. One notable 
exception is the 0{N) nonlinear sigma model in the limit of large A^, which is exactly soluble 
in an expanding = 1 FRW universe. In order to be reasonably self-contained, we shall 
follow Davis et al. |jTB[ closely in this section (note however the difference in normalizations 
of the random variables and the definition of the power-law exponent a). 

First, consider the lagrangian for the 0{N) fields 0", a = 1 ... A^: 

>c = ^v^rv^r-^(0)- (36) 

where we raise and lower indices with the background metric o?ri^y. If the 0{N) symmetry 
is broken to 0{N — 1), the potential will have a nonzero minimum, ^(0o) — 0; usual 
example is the broken "phi-fourth" potential, V{(f)) = A(0^ — 0o)^ with 0^ = With this 

sort of potential, the (p"' will have A^ massless modes corresponding to angular excitations, 
and one massive mode corresponding to radial excitations of the A^-vector 0. In the low- 
energy or strong-coupling limit, we can represent the potential of the massless modes with 
a lagrange multiplier term, 

£ = ^V.r'^'-r + - 0^) (37) 
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Thus, the fields ip"" are fixed to the (A^ — l)-sphere vacuum manifold: = 0q. The 

fields begin randomly distributed before the phase transition; as they come inside the Hubble 
radius, they come into causal contact with one another and order themselves to minimize 
the gradient energy (the first term of the Lagrangian) on those scales. 
The equation of motion is 

ro ro 

where the second equality follows from the vacuum manifold constraint enforced by the 
lagrange multiplier p2|] (and thus is only strictly true for the nonlinear sigma model, but 
not in models with a "physical" potential). At large A^, we can replace the bilinear V^^^V^c/)* 
with its (ensemble or spatial) average and the scaling ansatz 

{V,<P'W4>') = S^. (39) 

with a constant S. If this ansatz holds, the "density" of the stiff source is also proportional 
to and the dynamics are scale-invariant with respect to the horizon size (or Hubble 
radius). We shall see later that this choice is self-consistent. 

Before the phase transition (at temperatures T > 0o), the potential (or lagrange multi- 
plier) term in the lagrangian is irrelevant and we expect 0^ ~ due to thermal fluctations. 
After the phase transition, the field is pinned to the vacuum manifold. In the large- A^ limit, 
the distribution of the individual components 0" becomes a Gaussian peaked around (f)"' = 
with a variance given by (0^) = 0q (see Appendix). Therefore, we simply Fourier transform, 

0'^(x, r/) = J d'k a^0,(r^)e^'^- (40) 

where the prefactor enforces the constraint 0^ = 0o ^ind the are Gaussian random vari- 
ables of mean and normalized to unit variance for later ease of calculation. They are 
uncorrelated for unequal a and k; for example, 

(«)=5'^''5=^(k + q), (41) 

where 5"* and S^(k) are Kronecker and Dirac deltas, respectively, and higher (even) order 
correlations are sums of appropriate products of this two-point function; averages of odd 
numbers of the vanish. 

In a universe with scale factor a (x r]" and the scaling ansatz, Eq. (^) becomes the linear 
equation 

<I>1 + —<P'k + k'4>k = --,<Pk. (42) 

This has solution 

0,(ry) ^ k^'/'fikr^y, fix) = -L^i/2-j^^^(^) (43) 
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where the order of the Bessel function (and the value of the constant S) has been chosen 
so that the solution has a white-noise power spectrum on superhorizon scales, and A/q is 
a normalization constant. The presence of the simple form f{kri) indicates the expected 
scaling nature of the solution: kr] is the ratio of the horizon (r^ oc l/'H) to the length scale 
(A;~ 1/A). 

Given this solution, we can calculate the stress energy Q^^, required for the formalism of 
Section |T[ Most generally. 



e 



-9 



dg- 



(44) 



In this case, the second term vanishes and the "potential" (lagrange multiplier) term is zero 
when the vacuum manifold conditions are satisfied. So, 



Particular components are 



(45) 



(46) 
(47) 
(48) 



The solutions for the 0" can be inserted into these expressions to calculate the components 
in this model. We find that 



(Boo) oc — 



7]^ 



(49) 



in a matter-dominated universe, so the scaling solution holds. In a radiation-dominated uni- 
verse, there is also a factor of ln{r]/'r]i), so the scaling is modified by a logarithmic term. We 
also find that the "pressure," (0) vanishes in a matter-dominated universe and is (1/3) (0oo) 
during radiation-domination. We shall calculate Qoi below when considering density fluctu- 
ations. 



B. Density Perturbations and Power Spectra 

When the sigma-model solution, Eq. (^3]), and the fourier transform of the expression for 
the stress-energy from Eq. (|47|) are inserted into Eq. (|29|) above for the linear perturbation 
to the matter density, the result is 

5(k) = ^r^'J dr^^.Qo^ik) = ^r,' j ^/^gq-ka^^^ j dr^<P\^_^\{ri)Uv) (50) 

where we have assumed the power law exponent 7 = 0, or / drjOoi ^ const in Eq. (p9D, to 
be justified below. This form enables the computation of n-point power spectra. 
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In general, a function of the form {S(k)"') contains the product of 2n Gaussian random 
variables a^; thus even though the odd moments of the field may vanish, the corresponding 
moments of the density field may be nonzero, and thus may be non-Gaussian even to linear 
order. 

In addition, many of the terms in the expansion of the order n moment of the will 
vanish due to the form of the integrand above. For example, the power spectrum (5(k)5(k')) 
contains a term with 

k • qk' ■ q'5'^'^5"'"'5=^(k)5=^(k')A^"^ (51) 

where the Kronecker and Dirac deltas result from the mean of four of the from the 
previous equation. This term vanishes due to the presence of k'(5^(k') = 0. Similarly, in any 
moment, any term which contains a combination of Kronecker deltas that do not contract 
overall to will also contain a factor like k'5^(k') and will thus integrate to zero (this factor 
is effectively (5k) = 0); the remaining terms behave like 

k ■ qk' ■ q'5'*''*(5'^'"(5=^(k + q')53(k' + qjN'^ (52) 

in the calculation of the (2-point) power spectrum. In this case, the Kronecker deltas contract 
to one power of and the Dirac deltas are equivalent to 5^(k + k') after the integrals have 
been performed. For an n-point spectrum this works out to x A^~" = l/iV"~^, where 
the first factor comes from the contraction of the Kronecker deltas and A^"" from the 
n prefactors of 6(k) in Eq. (0). Moreover, when all the remaining integrals have been 
completed, there will still be a delta-function to enforce the requirement that the sum of the 
kj be zero as expected. Therefore, the moment will behave like 

(5(kO---5(k„))cx (^ry^y^5^(Ek.). (53) 

In fact, we can proceed further without actually doing any calculating. In addition to 
the factors in the previous equation, the nth moment will contain an integral over d^q of n 
quadratics in q and the kj as well as a product of n integrals of the form / dr] (f)p(j)p' , where 
p,p' are the lengths of linear combinations of q and the kj. As in Davis et al. [|16|, this 
integral can in general be written as a function of the configuration and overall scale {i.e., 
k = ki) of the kj by writing q = fcu and 1] = s/k, and all of the remaining kj = k{\<ii/k). 
Furthermore, we use the scaling of the (f) fields, 

UV) = p-'^'fiPV), <P'piv) = p-"^f\pv) (54) 

in the r] integral. We can then pull out all of the factors of scale to be left with /c^~" in front. 
Thus, the nth moment (factoring out the momentum-conserving 5 function) is given by 

P(ki, . . . , k„_i) = ( -^V^ 1 j;j—^9n{kr], configuration), (55) 

where "configuration" refers to the shape (but not the scale) of the n-sided polyhedron 
defined by the kj. (For n > 4, the irreducible moment P„ differs from the reducible (^i • • • 
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by terms that vanish unless two of the wavevectors satisfy kj + kj = 0; for continuous k 
this is a set of measure zero, and in any case does not apply to the configurations usually 
examined when considering data.) For suitably late times and subhorizon scales [ki] oo), 
the function gn should go to a constant, dependent only on the configuration of the figure. 

In particular we find that the subhorizon power spectrum is Harrison- Zel'dovich, P{k) oc 
k, the bispectrum -B(ki, k2) is independent of scale, and all higher moments go as a negative 
power of k] that is, they decrease with decreasing scale. However, this calculation is only 
valid for scales which have entered the horizon in the matter-dominated era. Thus, the 
spectrum of perturbations on scales smaller than the horizon at matter-radiation equality 
cannot be calculated from these formulae. 

Thus far, we have neglected the fact that we must take the large- limit of these quan- 
tities. If we normalize them to the (two-point) power spectrum P{k) through a quantity 
such as J3 (the volume integral of the two-point correlation function) or as, we get a finite 
numerical value for the combination G(j)1/\/N ^ 3 x 10~^cr8 |jl6[ (where (t| is the variance of 
the mass distribution in spheres of 8 /i~^Mpc; the variance in the number density of galaxies 
on that scale is measured to be cr8(gal) = 1 [^). The scaling of this value, (pl oc \fN is a 
consequence of the analytic calculations reproduced in this section for large N . However, 
the numerical value comes from various simulations. Comparison of the normalization (to 
(Tg) for specific values of < 6 in simulations (which explicitly include topological-defect 
field configurations not present for large A^), indicates that this value for 0o may be as much 
as a factor of 10 too high. Although the calculations reproduced in this section show that 
we should expect 0o °^ in the large- iV limit, for these low values of that have been 
numerically investigated, such scaling has not yet been reached 1^^. In any case, none of 
these values are expected to be known to better than about 50% fl^ (The low-A^ cases are 
in the process of being solved by using exact Green's functions rather than direct integration 
of the equations of motion presented here. This is expected to give more accurate results 



In the large- A^ limit, these quantities decrease with A^ for all > 3. However, the leading 
term in 1/A^ causes a non-Gaussian distribution to this order, and we expect these calcula- 
tions to be at least approximately valid for A^ greater than a few, since the character of the 
sigma-field is approximately the same in the absence of topological defects such as strings 
(A^ = 2) or textures (A^ = 4). 

The asymptotic value of (72 for large k (small scales) is given by 



Thus, the n-point spectra work out to be 




(56) 




(57) 



as in where 




(58) 
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We have already assumed that fc?] S> 1 so we may integrate by parts and ehminate surface 
terms to consohdate the various integrals /(a, b; kt]) = I{a, h) that appear (see below and 
Fig. (|1|)). The corresponding value for the bispectrum, g^,, is 



^3(v) = j rf3uu-k/(|k-u|,M) 



[1 + 2v ■ k - 2u ■ k + t;^ - 2u • v)(2u • V - 2v • k - w^) X 
/(|k + V - u|, |u - k|)/(M, |k + V - u|) + 



[v"^ + 2u ■ v)(2u ■ k + 2u ■ - l)/(n, |u + v|)/(|u - k|, |v + u|) 



(59) 



where v = \s~2/ki embodies the dependence upon the configuration of the wavevectors. For 
an equilateral triangle, k2 = ki = k and u ■ v = — 1/2. 

Note that these quantities do have a residual scale dependence in the form of the upper 
limit of the innermost integrals I. In Fig. (|l]), we plot numerical integrations of g2 and 
gs as a function of k where we integrate from horizon crossing [kr] = 1) to the present 
epoch [kr] = krjo) at each wavenumber. Both of these quantities reach their asymptotic 
values {g2 ~ 12, (73 ~ 1.6) by k~^ ^ 10^ h~^Mpc, well outside of the range of our ability to 
reliably measure higher-order correlation functions. This quick approach to constant values 
indicates that, as suspected, / drjdiQoi ~ const, so we are justified in assuming the power- 
law exponent 7 in Eqs. (^) and (^2]) above — the situation mimics a simple growing 
mode. Physically, this is simple to understand: as a given mode enters the horizon, it comes 
into causal contact fairly quickly, on the order of a Hubble time. After that, the field is 
roughly static, so its momentum density Bqi is small, and the matter perturbations are 
only generated on scales near the horizon. Subsequently, they evolve under the influence of 
small-scale Newtonian gravity alone. (However, it is not hard to imagine some other kind 
of stiff source with "ordering physics" that continues to be active on smaller scales which 
would result in density fluctuations that do not mimic the Newtonian growing- mode result.) 

One traditional way to characterize the bispectrum is through the ratio 

5(ki,k2) , , 

Q = ^ ^' ^' (60) 

^ P,P2 + PlP3 + P2P3 ^ ' 

which for an equilateral triangle becomes Q = B[k)/3P{k)'^. The usual hierarchical result 
from the Newtonian analysis with Gaussian initial conditions is Q = const from second-order 
perturbation theory. With the non-Gaussian initial conditions from the sigma model, we 
instead have Q oc from linear theory. 

The question remains, then, on what scales do we expect these linearly evolved mo- 
ments to dominate over higher-order effects? From Section |II E| above, the leading nonlinear 
contribution to the higher moments is 



/ \2{n-l) 

PnAk) = QnPr' = Qn f ^^^^ j (^^^^O""', (61) 



where this equation is actually only correct for n = 3 with equilateral triangles, where 
Qa = SQ. (For n > 3, the right-hand side of this equation should be a sum over the 
different "tree-level" graphs connecting the labelings of the wavevectors, each term with a 
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corresponding constant Qn-a |18|, [T0|J26[] . However, for analyzing the scale-dependencies here, 
the form above shall suffice.) 

Consider the quantities dn = Pn/P'z^'^ which normalize the ra-point spectra to |(5(k)|". 
We see that the linear contribution is, of course, constant, 

duM = ^A:3(i-"/2)iVi-"/2 (62) 

92 

whereas the nonlinear contribution grows with time. 



dn,ni = Qn ( -^g2kA ^ (63) 



At a given wavenumber, the nonlinear contribution grows with time with respect to the linear 



contribution. A similar result for spatial correlation functions was derived for n = 3 in pT 
where they also included the effects of a primordial four-point function (whose contribution 
to the skewness also grows with time). 

The ratio of the nonlinear contribution to that of the linear evolution is 



n,nl ^ 



Pn. . 1 i 



n,lin 



r^A2 -|"-2 n-1 



which goes as a positive power of k for n > 2, so the nonlinear terms become more important 
on small scales, as expected. This is a function of krj, the ratio of the length involved to the 
horizon size, so the scaling nature of the solution is preserved. 

The scale of this crossover from the domination of the linear to the nonlinear evolution 
of the fields depends on the values of the Qn for a given configuration of wavevectors and 

(although the powers of A^ explicitly cancel in the above expression, recall that the 
normalization of 00 is such that factor in brackets above is constant). Crossover thus occurs 
at scale 

«„M N^'i^) . (65) 



Actually, the full evolution of the three-point function depends on the four-point function 
(or trispectrum T(ki, k2, ks)) as well which we do not explicitly calculate here, but we 
expect this crossover scale to be independent of it, since from above we have T oc k~^ so 
the trispectrum falls on small scales, at least for those greater than that of matter-radiation 
equality. Moreover, the value of n only enters these equations through the prefactor involving 
the Qn and Qn-, which are not expected to vary greatly with n. So, for increasing n, the 
n-dependent part of the first factor in parentheses will be come less important, and the 
crossover scale becomes approximately 

k-^ ^ mT^'" (y^) ; (large n). (66) 



If we naively plug the asymptotic values of g2 ~ 12 and (73 ~ 1.6 into Eq. (|65D, with the 
hierarchical value for equilateral triangles of = 3(5 = 34/7 we find 
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fc-j^ ^ 219N^/^al^\~^Mpc; {n = 3) (67) 

for G(f)l/\/N ^ 3 X lO^^cTs- If we let this quantity be a factor of 10 smaller as indicated 
above, then the crossover scale is reduced to k~^^ ^ 69N^^'^al^'^ h~^Mpc. To calculate the 
scale at which the primordial four-point function or trispectrum becomes significant, we will 
use the hierarchical pattern T ^ 16-RP^ for tetrahedral configurations, with R ~ Q^, [|| so 
Q4 ~ IQQ^ ~ 5 (although observations [^] are typically R ^ 1.) Thus, 

k-^^ ^ (17 - 54)^4"^/^ A^^/Vg/^/i-^Mpc; {n = 4) (68) 

spanning the range of possible normalizations of (f)Q, where we have left the integral 
unevaluated, but expect it to be of order unity. 

Even this scale, unfortunately, is just beyond the largest considered in surveys for higher- 
order power spectra, which have gone out to k~^ ^ 20 /i^^Mpc ||2^. We expect only the 



non/mear contribution to the three-point function to be important on the smallest observable 
scales. If instead we use the second, large-n form of k^i, we find 

fe'i^ ~ (11 - 3Q)N^/^al^^h-^Mpc; (large n) (69) 

This is more possibly in the range of observed and observable scales, but still (for ^ 6 
where we expect to believe these results) barely on the edge of current observations of the 
bispectrum and trispectrum. 

Thus, the usual Newtonian analysis of higher-order correlation should suffice there. Un- 
fortunately, that means that it will be more difficult to distinguish the nonlinear sigma model 
from simple Gaussian theories on small scales using the mass distribution alone. However, 



as discussed in ^7j, the large-scale normalization of the nonlinear sigma model to the COBE 
results may require an unreasonably large bias in order to match the galaxy distribution. 

More problematic is the calculation on scales that entered the horizon during the 
radiation-dominated period. A complete calculation requires the numerical integration of 
the equations of motion of the sigma model fields and the radiation and matter fluids as in 
T7| . They present a possible parametrization of the transfer function which takes the power 
spectrum from P{k) (x k on large scales to its small-scale form (usually P{k) oc k~^, pos- 
sibly modulated by the logarithmic growth of modes in the radiation-dominated universe). 
However, the equivalent form has not been calculated for the three-point function. However, 
we do not expect the linear evolution to dominate on small scales. Higher order or fully 
nonlinear effects tend to dominate on the smallest scales. In order for the linear correlations 
to be observable below the scale of matter-radiation equality, the linear bispectrum would 
have to increase relative to the nonlinear contribution on these scales. Moreover, because of 
the form of Eqs. ( P3a| ), the relative importance of higher order effects with respect to linear 
effects is again expected to be the same as in the purely Newtonian case. 



C. Comparison with Observations 

As we have seen, the nonlinear sigma model is not expected to give non-Gaussian results 
for higher-order correlation functions on as-yet observed scales. So far, all observations are 
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consistent with Newtonian hierarchical results with Gaussian initial conditions. Baumgart 
and Fry [^] actually analyze the n-point power spectra for two galaxy surveys and find 
that Q{k) = B/3P^ ~ 4/7 (for equilateral triangles) for ^ 10 /i~^Mpc. This is exactly 
the value expected from the Newtonian theory with Gaussian initial conditions. They also 
calculate R{k) = T/IQP'^ ^ 1 (for regular tetrahedral configurations), a result which is 
weakly consistent with the hierarchical expectation. 

More recent results of analyses of the IRAS dataset are also consistent with the hi- 



erarchical expectations. Bouchet et al. [27] have calculated the volume-averaged n-point 



functions for n <b and again find that the hierarchical form holds: oc ^2 ^- Moreover, 



Nusser et al. [^] have reconstructed the initial one-point probability density function for 
the smoothed IRAS density field (using an algorithm based on Newtonian gravity and the 
Zel'dovich approximation) and found it to be consistent with a Gaussian distribution on 
scales smaller than about 70 /;,~^Mpc. 

It may be possible to examine the n-point functions at still larger scales using two- 
dimensional angular correlations . However, an angular analysis introduces errors due to 
uncertainties in the luminosity function of galaxies and consequently the survey's selection 
function, especially on the large scales we are most interested in. Nonetheless, Gaztafiaga 
pUI has estimated the n-point angular correlation functions for the APM dataset and found 
that the results are consistent with an initially Gaussian distribution for scales out to about 
30/i^^Mpc; however, this is still below the scales required to see any possible effects of a 
sigma-model field. 

Precise measurements of the bispectrum and trispectrum out to scales of /c^^ ~ 
100 /i~^Mpc should be sufficient to begin to more explicitly test this model. We should 
expect the next generation of survey results for higher-order correlations to at least begin to 
probe the scales on which primordial non-Gaussianity from a sigma-model source should be 
observable, at least for moderate values of A^. That we do not yet observe a marked increase 
in the three-, four-, or five-point functions (or their spectral counterparts) on the largest 
scales implies that the density field on those scales remains dominated by Newtonian evolu- 
tion. If this observational trend continues, then any crossover to linear evolution will be on 
still larger scales, and we will be able to make stronger statements about the normalization 
or existence of any sigma-model field {e.g., using Eq. (|65D). 

We re-emphasize the fact that observations of the scaling hierarchy of correlation function 
do not imply Gaussian initial conditions, even on the scale of the observations. As we have 
seen here, the correlation functions on small scales are dominated by the contributions from 
the initially Gaussian component of the density field, which evolves nonlinearly into the 
quasi- Gaussian scaling hierarchy. This behavior is expected to be generic to any scenario 
with non-Gaussian initial conditions, as the nonlinear evolution should in any case dominate 
on small scales. 

Other observational ramifications of the sigma-model may be more damning even to- 



day. As Pen et al. [|l^] point out, normalization of the microwave distortions produced in 



these models to the COBE observations require what may be an uncomfortably high bias, 
h = l/o"8 ~ 2/h. Moreover, the detailed fit of the power spectrum of mass fiuctuations 
(determined from numerical simulations) to the measured QDOT galaxy distribution seems 
to require a bias of order 6 on scales of 20 h~^Mpc. 
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IV. CONCLUSION 



We have developed a formalism to examine the quasi-nonlinear behavior of perturba- 
tions in a universe with a "stiff" source. In particular, we have shown that their evolution 
is formally very similar to those encountered in a purely Newtonian analysis with an ini- 
tial perturbation spectrum, but without a source. In the particular case of the nonlinear 
sigma model, we find that the behavior of higher correlation functions of the density field 
on extreme subhorizon scales is exactly the same as that which occurs with primordial adi- 
abatic perturbations. On larger, currently unobservable, scales, two effects occur that might 
differentiate the sigma model. First, the sigma model (at least at moderate A^), has an 
initially non-Gaussian distribution, which might be observed on scales of several hundred 
megaparsecs. Furthermore, we find that the density perturbations are actually being created 
on scales comparable to the horizon, where the self-ordering physics of the source is occur- 
ring. Unfortunately, the prospects of observing effects on these scales are slim. However, 
it may be possible to invent models in which the equivalent self-ordering occurs on smaller 
scales as well. Unfortunately, the chief candidates — defect theories like cosmic strings — are 
ill-suited for this analysis, because the gravitational effects of defects generally occur on too 
small a scale. Therefore, a quasi-nonlinear analysis is insufficient and full-blown numerical 
simulations must be performed. 

There are several possible extensions to the work we have presented here. We have 
concentrated on scalar perturbations, because those result in the readily-observable density 
fluctuations, even when higher-order terms in the matter variables are present. By consider- 
ing vector and tensor perturbations, we can go beyond this and calculate quantities like the 
full velocity field (including any possible vortical part), as well as the tensor perturbations 
(gravitational radiation) from the stiff source. 

If we consider a universe filled with both matter and radiation fluids, it should be pos- 
sible to use this formalism to make accurate calculations of microwave anisotropies in these 
theories. If we have also calculated the gravitational radiation, then we can separate out 



the scalar and tensor components [31 



Finally, it should be possible to extend the cosmological post-Newtonian suite of ap- 
proximations pO| to include the possibility of a stiff source. Because this involves going to 
a higher order in the metric perturbation, things become considerable more complicated, 
but such approximations are often useful even when perturbations have become extremely 
nonlinear. 
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APPENDIX: DISTRIBUTION OF THE FIELDS IN THE 0{N) MODEL 

In the nonlinear 0{N) sigma model, the fields 0" are constrained to lie upon the vacuum 
manifold, an {N — l)-sphere of radius 0o, where they will be distributed uniformly in the 
absence of any further symmetry breaking. Thus, the probablity density will be proportional 
to the solid angle on the (A^ — l)-sphcrc. 

We choose A'" — 1 polar coordinates {ip, 9i, . . ., 9n^2) such that the solid angle is 

dQjsf = dip sin 9i d9i--- sin' 9i d9i - ■■ sin^~^ 9^-2 d9N-2- (Al) 

Thus, the distribution of fields in polar coordinates is given by 

p{(p, {9i}) dip d9i - ■■ d9M-2 oc d9,N- (A2) 

We are concerned, however, with the distribution of the 0", the Cartesian components of 
the fields. In this coordinate system, there is always one component given by ^ = 0"/0o = 
cos^Ar-2, and, due to the symmetry of the system, we can calculate the distribution of any 
single Cartesian component: 

p{z) = j p{ip, {9i}) dip d9i - ■ ■ d9N-2 5{z - cos 9n-2) 
— ^' j dQ:N5{z — cos 

= A' J dip J d9i sin ^1 • • • y" sin^-^ ^^.3 J d9N-2 sin^-=^ 5{z - cos 9^-2) 

= 2t:A! j d9i sin ^1 • • • y d9N-3 sin""-^ 9n-3 f dy (1 - y2)(^-3)/25(^ _ ^) 

= A{1 - ^2)(^-3)/2 (A3) 

where A, A' are constants determined by the requirement that / dzp{z) — 1, and we have 
made the change of variables y = cos^jv-2- This gives 



/ X , 1 T(N/2) ( o\(Af-3)/2 , 



This distribution has mean (^) = and variance {z^) = 1/N (which can be calculated 
directly or seen from the constraint J2<P"'<P"' — (Po)- The higher moments are given by 

_ 1 r(iV/2)r(m/2 + l/2) 
^ " A r(iV/2 + m/2) ^^^""'^ ^^^^ 

which behave as 0(1/A^™/^) for large N; the odd-m moments vanish due to the symmetry 
of the distribution. For large A^, 

P(z)^ll(l-zf'\ (A6) 
We wish to compare this to a Gaussian distribution with some width a, 
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1 




r(m/2 



+ 



1/2) (m even). 



(A7) 



For small and large A^, 




V27r(T2 



1 




(A8) 



which implies = to first order in as expected from {z^) = 1/N and the 

expression for the higher moments of the Gaussian distribution at large A^. For z ^ 1, both 
p{z),g{z) — > for large A^ or small a (although the ratio p{z)/g{z) can be quite large for a 
finite value of A^; this is because the Gaussian is normalized over the interval (— oo, +cxd)). 
Moreover, in the limit of A^ — oo, both p{z) and g{z) approach the Dirac 5-function. In 
Figure || we show p{z) and g{z) for various values of A^. Note that for A^ < 5, the departures 
from Gaussianity are significant. 

Thus, the distribution of the Cartesian components of the field approaches that of a 
Gaussian with o"^ = 1/A^. For large A^, then, the components of the field become more 
and more strongly peaked about z = 0. Note that this distribution is applicable only on 
super horizon scales, where the ordering dynamics of the field are unimportant. Within the 
horizon, the fields will continually organize in order to minimize their gradient energy as 
new scales enter the causally-connected region. 
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FIGURES 



FIG. 1. The integrals g2{k) (solid) and g^ik) (dashed) which are the prefactors in the cal- 
cuation of the linear contribution to the power spectrum P = P2 (x and bispectrum (for 
equilateral triangles) -B = P3 oc g^k^ respectively (Eq. (|55|)). We have integrated from hori- 
zon-crossing [kr] = 1) to the present day [kr] = krjo). In each case the asymptotic value is reached 
by A;"^ ~ 1000/i~^Mpc. Note that we have already assumed krj ^ 1 in writing down the expres- 
sions for the gn, so the details of the approach to the asymptotic values (32 — > 12.2 and 53 1.6) 
should only be taken as indicative. 

FIG. 2. A comparison of ^(2: = /(po) (solid), the actual distribution of fields <j)°' in the 0{N) 
sigma model, with g{z) (dashed), the corresponding Gaussian distribution of mean and variance 
for N = 3,5,10 and 50. 
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